library(flowCore)
library(flowStats)
library(flowViz)
library(ggplot2) 
##library(lessR)
##install.packages("lessR")
##library(gridExtra)
##library(multicore)




if(.Platform$OS.type=="windows"){
  root<-"D:\\mydocs\\2013\\dev\\ProteinSW\\"
  data.dir<-"D:\\mydocs\\2013\\dev\\ProteinSW\\"
}else{
  root<-"/home/haseong/alpha/dev/FACS/"
  data.dir<-"/home/haseong/alpha/dev/FACS/ProteinSW/"  
  data.dir2<-"/home/haseong/alpha/dev/FACS/ProteinSW/20130927/"  
}
setwd(root)


tick_locations = c(0, 256, 512, 768, 1024)  
tick_labels = 10^(tick_locations/256) 



prepro.facs<-function(x, ncells=10000, view.plot=T){
  z<-list() ## for saving results
  xx<-exprs(x)  
  ## delete false
  del.idx<-which(xx[,"SSC-H"]<=50)
  del.idx<-c(del.idx, which(xx[,"FL1-H"]<=0))
  xx<-xx[-del.idx,]
  
  ## select 10000 cells from
  ##ncells<-15000
  ## select the most dense area
  fsc.density<-density(xx[,"FSC-H"])
  fsc.max<-fsc.density$x[which.max(fsc.density$y)]
  ##fsc.max<-256
  ssc.density<-density(xx[,"SSC-H"])
  ssc.max<-ssc.density$x[which.max(ssc.density$y)]
  ##ssc.max<-356
  d<-sqrt((xx[,"FSC-H"]-fsc.max)^2+(xx[,"SSC-H"]-ssc.max)^2)
  o<-order(d)
  sel.idx<-o[1:ncells]
  sel.xx<-xx[sel.idx,]
  cat(fsc.max, ssc.max, "mean:", mean(sel.xx[,"FL1-H"]), "dim:", dim(sel.xx), "\n");flush.console()
  
  z$fsc.max<-fsc.max
  z$ssc.max<-ssc.max
  ##z$fsc.mean<-mean(sel.xx[,"FSC-H"])
  ##z$ssc.mean<-mean(sel.xx[,"SSC-H"])
  
  ### ----- plot ------
  if(view.plot==T){
    par(mfrow=c(3,2))
    ## FSC vs cell density
    plot(density(xx[,"FSC-H"]), col="#000000", ylab="density", xlab="FSC", xaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    
    ## SSC vs. cell density
    plot(density(xx[,"SSC-H"]), col="#000000", ylab="density", xlab="SSC", xaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    
    ## FSC vs. SSC and the automatically selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"SSC-H"], pch=20, col="#00000001", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xaxt = "n", yaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    points(sel.xx[,"FSC-H"], sel.xx[,"SSC-H"], pch=19, col="#0000ff03", xlab="FSC", ylab="SSC", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    ### FL1 of the selected cells vs. their density
    plot(density(xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", bty="n", main="", xaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    
    ## FSC vs. FL1 without the selected 10000 cells
    plot(xx[,"FSC-H"], xx[,"FL1-H"], pch=20, col="#555555", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. FL1)", xaxt = "n", yaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    axis(2, at = tick_locations, labels = tick_labels)
    points(sel.xx[,"FSC-H"], sel.xx[,"FL1-H"], pch=19, col="#0000ff", xlab="FSC", ylab="FL1", main="Scatter plot (FSC vs. SSC)", xlim=c(0,500), ylim=c(0,800))
    
    ### FL1 vs. cell density
    plot(density(sel.xx[,"FL1-H"]), col="#000000", ylab="density", xlab="FL1", bty="n", main="", xaxt="n")
    axis(1, at = tick_locations, labels = tick_labels)
    #text(0, mean(density(sel.xx[,"FL1-H"])$y), paste(data.info[1, "effector"], data.info[1, "effector.level"], sep=" "), pos=4)
    
  }
  
  z$sel.xx<-sel.xx
  return(z)
  
}




x.neg1<-read.FCS(paste(data.dir, "pucbb-laci-neg1.001", sep=""), transformation=FALSE)  
x.neg2<-read.FCS(paste(data.dir, "pucbb-laci-neg2.002", sep=""), transformation=FALSE) 
x.neg3<-read.FCS(paste(data.dir, "pUCBBLacI.fcs", sep=""), transformation=FALSE) 
x.pos1<-read.FCS(paste(data.dir, "pucbb-laci-pos2.003", sep=""), transformation=FALSE) 
x.lib2<-read.FCS(paste(data.dir, "LacI2ndlibrary.fcs", sep=""), transformation=FALSE)  
##x<-read.FCS(paste(data.dir, "eGFP2.fcs", sep=""), transformation=FALSE)  
x.neg4<-read.FCS(paste(data.dir2, "lacNeg4sort.fcs", sep=""), transformation=FALSE)  
x.neg5<-read.FCS(paste(data.dir2, "laciNeg.fcs", sep=""), transformation=FALSE)  
x.lib<-read.FCS(paste(data.dir2, "lacLib.fcs", sep=""), transformation=FALSE)  


neg1<-facs.summary.plot(x.neg1)
neg2<-facs.summary.plot(x.neg2)
neg3<-facs.summary.plot(x.neg3)
neg4<-facs.summary.plot(x.neg4)
neg5<-facs.summary.plot(x.neg5)

lib<-facs.summary.plot(x.lib)
lib2<-facs.summary.plot(x.lib2)

pos1<-prepro.facs(x.pos1)







